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Abstract. Algorithms for the symboHc computation of conserved densities, 
fiuxes, generalized symmetries, and recursion operators for systems of nonlin- 
ear differential-difference equations are presented. In the algorithms we use 
discrete versions of the Frechet and variational derivatives, as well as discrete 
Euler and homotopy operators. 

The algorithms are illustrated for prototypical nonlinear lattices, includ- 
ing the Kac-van Moerbeke (Volterra) and Toda lattices. Results are shown for 
the modified Volterra and Ablowitz-Ladik lattices. 



1. Introduction 

The study of complete integrability of nonlinear differential-difference equa- 
tions (DDEs) largely parallels that of PDEs AC91, ASYOO, W98 . Indeed, as 
in the continuous case, the existence of sufficiently many conserved densities and 
generalized symmetries is a predictor for complete integrability. Based on the first 
few densities and symmetries a recursion operator (which maps symmetries to sym- 
metries) can be constructed. The existence of a recursion operator, which allows 
one to generate an infinite sequence of symmetries, confirms complete integrability. 

There is a vast body of work on the complete integrability of DDEs. Consult 
e.g. |ASYOO, HH03; for various approaches and references. In this article we 
describe algorithms to symbolically compute polynomial conservation laws, fluxes, 
generalized symmetries, and recursion operators for DDEs. The design of these 
algorithms heavily relies on related work for PDEs |GH97al IHG99L IHGCM98j 
and work by Oevel et al |OZF93| . 

There exists a close analogy between the continuous and discrete (in space) 
cases. As shown in |HM02L IMH02| . this analogy can be completely formalized 
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and both theories can be formulated in terms of complexes. The same applies for the 
formulation in terms of Lie algebra complexes in jW98j . This allows one to translate 
by analogy the existing algorithms immediately (although complications arise when 
there is explicit dependence on the space variable in the discrete case) . One of the 
more useful tools following from the abstract theory is the homotopy operator, 
which is based on scaling vectorfields, and goes back to Poincare in the continuous 
case. This operator allows one to directly integrate differential forms and can 
be straightforwardly implemented in computer algebra packages, since it reduces 
to integration over one scaling parameter. The discrete analogue as described in 
|KM02, MH02 does the corresponding job in the discrete case and we use it to 
compute fluxes. In this paper we do not explicitly use the abstract framework, yet, 
it has been a motivating force for the development of our algorithms. 

The algorithms in this paper can be implemented in any computer algebra 
system. Our Mathematica package InvariantsSymmetries.m GH97b computes 
densities and symmetries, and therefore aids in automated testing of complete in- 
tegrability of semi-discrete lattices. Mathematica code to automatically compute 
recursion operators is still under development. 

The paper is organized as follows. In Section [21 we give key definitions and 
introduce the Kac-van Moerbeke (KvM) iKvM75 , Toda |T8T| and Ablowitz-Ladik 
(AL) lattices |AL75j . which will be used as prototypical examples throughout the 
paper. The discrete higher Euler operators (variational derivatives) and the discrete 
homotopy operator are introduced in Section |3 These operators are applied in the 
construction of densities and fluxes in Section The algorithm for higher-order 
symmetries is outlined in Sectional The algorithms for scalar and matrix recursion 
operators are given in Section [7| and Section |51 The paper concludes with two more 
examples in Section |5| 

2. Key Definitions and Prototypical Examples 
Definition 2.1. A nonlinear (autonomous) DDE is an equation of the form 

(2.1) u„ = F(..., u„_i, u„, u„+i, ...), 

where u„ and F are vector-valued functions with m components. The integer n 
corresponds to discretization in space; the dot denotes one derivative with respect 
to the continuous time variable {t). 

For simplicity, we will denote the components of u„ by (u„, w„, m„, • • • ). For 
brevity, we write F(u„), although F typically depends on u„ and a flnite number 
of its forward and backward shifts. We assume that F is polynomial with constant 
coefficients. If present, parameters are denoted by lower-case Greek letters. No 
restrictions are imposed on the level of the forward or backward shifts or the degree 
of nonlinearity in F. 

Example 2.2. The AL lattice |AL75j . 

Un = [Un+l - 2u„ -I- Un-l) + U„W„(u„+l + U„-l), 

(2.2) Vn = -{Vn+l-'2.Vn+Vn-l)-UnVn{Vn+l+Vn-l), 

is a completely integrable discretization of the nonlinear Schrodinger equation. 

Definition 2.3. A DDE is said to be dilation invariant if it is invariant under 
a scaling (dilation) symmetry. 
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Example 2.4. The KvM lattice |KvM75j . 

(2.3) iln = UniUn+l - Un-l), 

is invariant under (t,u„) — > {X^^t, Au„), where A is an arbitrary scaUng parameter. 
Example 2.5. The Toda lattice |T8T] in polynomial form |GH98| . 

(2.4) W„ = U„_i - «„, Vri = v„{u„ - U„+i), 

is invariant under the scaling symmetry 

(2.5) {t, u„, Vn) (A"4, A-u„, A^w„). 

Thus, Un and i>„ correspond to one, respectively two derivatives with respect to t, 

dt dt'' 



(2.6) ' - -,::2- 



Definition 2.6. We define the weight, w, of a variable as the number of t- 
derivatives the variable corresponds to. 

Since t is replaced by j, we set w{-^) — 1. In view of H2.6|) . we have w(u,i) — 1, 
and w(wri) = 2 for the Toda lattice. Weights of dependent variables are nonnegative, 
rational, and independent of n. 

Definition 2.7. The rank of a monomial is defined as the total weight of the 
monomial. An expression is uniform in rank if all its terms have the same rank. 

Example 2.8. In the first equation of 12.4|l . all the monomials have rank 2; in 
the second equation all the monomials have rank 3. Conversely, requiring uniformity 
in rank for each equation in (|2.4I) allows one to compute the weights of the dependent 
variables with simple linear algebra. Indeed, 

(2.7) w{Un) + l = w{Vn), w{Vn) + I ^ w{Un) + w{Vn) , 

yields 

(2.8) w{u,,) = 1, w(w„) = 2, 
which is consistent with (|2.5() . 

Dilation symmetries, which are special Lie-point symmetries, are common to 
many lattice equations. Lattices that do not admit a dilation symmetry can be 
made scaling invariant by extending the set of dependent variables using auxiliary 
parameters with scaling. 

Example 2.9. The AL lattice is not dilation invariant. Introducing an auxiliary 
parameter a, hence replacing H2.2|l by 

iln = a(un+i - 2u„ + ?i„_i) + it„w„(u„+i + U„_i), 

(2.9) in = -a{Vn+l - 2Vn+Vn-l) - UnVniVn+1 + Vn-l), 

and requiring uniformity in rank, gives 

w(m„) + 1 = w{a) + 'w{Un) ^ 2w{Un) + w{Vn), 

(2.10) w{Vn) + l = w{a) + w(Vn) ^ 2w{Vn) + w{Un). 

Obviously, 

(2.11) w{u„) + w{vn) ^ w{a) ^ 1. 
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Several scales are possible. The choice w{un) = w{vn) = ^,w{a) = 1, corresponds 
to the scaling symmetry 

(2.12) {t,Un,Vn,a) (A"4, A^u„, A'w„, Aa). 

Definition 2.10. A scalar function p„(u„) is a conserved density of (|2.1() if 
there exists a scalar function J„(u„) called the associated flux, such that 

(2.13) Dtp„ + A J„ = 
is satisfied on the solutions of (|2.1|l |093| . 

In H2.13|l . we used the (forward) difference operator, A = D — I, defined by 

(2.14) AJ„ = (D-I)J„ = J„+i- J„, 

where D denotes the up-shift (forward or right-shift) operator, DJ„ = Jn+i, and 
I is the identity operator. Operator A takes the role of a spatial derivative on 
the shifted variables as many DDEs arise from discretization of a PDE in (1 + 1) 
variables. Most, but not all, densities are polynomial in u„. 

Example 2.11. The first three conservation laws for H2.3|l are 

(2.15) Dt(ln(u„)) + {Un+l + Un) - [Un + Un-l) = 0, 

(2.16) Dt(u„) + (u„+iu„) - (m„u„_i) = 0, 

(2.17) T)t{]^ui + U„U„+i) + U„W„+i(u„+i + Un+2) - Un-lUn{Un + U„+l) = 0. 

Example 2.12. For the Toda lattice (|2.4|l the first four density- flux pairs are 



(2.18) 


p(0) 


= ln(z;n), 






(2.19) 


p(i) 

rn 






= Vn-1, 


(2.20) 


rn 


1 2 




^n^n— 1 1 


(2.21) 


p(3) 
rn 


1 3 . 


7(3) 





The above densities are uniform of ranks through 3. The corresponding fluxes 
are also uniform in rank with ranks 1 through 4. In general, if in 12.13|l rank p„ = R 
then rank Jn — R -\- 1, since w{T)t) = 1- 

This comes as no surprise since the conservation law ()2.13|l holds on solutions 
of (|2.1|l . hence it 'inherits' the dilation symmetry of (|2.1() . 

In Section 13 we will give an algorithm to compute polynomial conserved densi- 
ties and fluxes and use H2.4|l to illustrate the steps. Non-polynomial densities (which 
are easy to find by hand) can be computed with the method given in HH(^. 

Definition 2.13. Compositions of D and D^^ define an equivalence relation 
(=) on monomials in the components of u„. All shifted monomials are equivalent. 

Example 2.14. For example, Un-iVn+i = Un+2Vn+4 = Un-sVn-i- 

Definition 2.15. The main representative of an equivalence class is the mono- 
mial of the class with n as lowest label on any component of u (or i; if m is missing). 

Example 2.16. For example, UnUn+2 is the main representative of the class 
{• ■ • , u„_2W„, u„_iM„+i, u„M„+2, Un+iUn+3, ■■ ■}■ Usc Icxicographical ordering to 
resolve confiicts. For example, m„i'„+2 (not Un-2Vn) is the main representative 

of the class {• ■ • ,U„_3l)„_i, • ■ • ,U„+2t'n+4, • ■ ■ }• 
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Definition 2.17. A vector function G(u„) is called a generalized symmetry 
of (|2.1() if the infinitesimal transformation u„ — > u„ + eG leaves H2.1|l invariant up 
to order e. Consequently, G must satisfy j093j 

(2.22) DtG = F'(u„)[G], 

on solutions of (|2.1|l . F'(u„)[G] is the Frechet derivative of F in the direction ofG. 
For the scalar case {N = 1), the Frechet derivative in the direction of G is 

(2.23) F'{un)[G] = ^F{ur, + eG)U=o = V tt^D^G, 

k 

which defines the Frechet derivative operator 

(2.24) F'[u^)^Y.T^ 

In the vector case with two components, u„ and u„, the Frechet derivative operator 
is 



(2.25) F'(u„) 



I l^k du^+k^ l^k dv^+f,^ 
_dF2_Y\k _dF2_Y\k 



Applied to G = (Gi G2)^ , where T is transpose, one obtains 

(2.26) ^^/(u„)[G] = V-^D''^Gi+V^D^-G2, z-1,2. 

In (|2.23() and 12.2611 summation is over all positive and negative shifts (including 
fc = 0). For fc>0, D''=DoDo---oD(fc times). Similarly, for fc < the down-shift 
operator is applied repeatedly. The generalization of (|2.25(l to a system with 
N components should be obvious. 

Example 2.18. The first two symmetries of H2.3|l are 

(2.27) G(i) = «„(«„+! -u„_i), 

(2.28) G*^^' = M„M„+i(m„ + M„+i + Un+2) - Un-lUniUn^2 + + Un) ■ 

These symmetries are uniform in rank (rank 2 and 3, respectively). The sym- 
metries of ranks and 1 are both zero. 

Example 2.19. The first two non-trivial symmetries of 1)2.4(1 are 

(2.29) G(i) = 



Vn-l - Vn 
Vn{Un - Un+l) 

and 

4+1 -ul + Vn+i -W„-l) 



(2.30) G(2) 



(2) 

The above symmetries are uniform in rank. For example, rankG-j^ = 3 and 
(2) 

rankGa = 4. Symmetries of lower ranks are trivial. 

An algorithm to compute generalized symmetries will be outlined in Section |^ 
and appHed to H2.4|l . 
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Definition 2.20. A recursion operator TZ connects symmetries 

(2.31) G(J+'') =7^G(^■^ 

where j = 1,2,..., and s is the seed. In most cases the symmetries are consecutively 
Hnked (s — 1). For A^-component systems, 7^ is an iV x A matrix operator. 
The defining equation for 7^ [0931 IW98j is 

dTZ 

(2.32) Dt7^ + [TZ, F'(u„)] = ^ + 7^'[F] + 7^ o F'(u„) - r'(u„) o 7^ = 0, 

ot 

where [ , ] denotes the commutator and o the composition of operators. The 
operator F'(u„) was defined in l|2.25|l . 7^'[F] is the Frechet derivative of TZ in the 
direction of F. For the scalar case, the operator TZ is often of the form 

(2.33) TZ = U{un) OiiB - I)-\D-\l,D)V{un), 
and then 

(2.34) 7^'[F] =^(D'=F)-^^Ol^ + ^ {/0(D''^F) 



k k 

For the vector case, the elements of the N x N operator matrix TZ are often of the 
form TZij = J7ij(u„) ((D — I)~^, D~^, I, D) (u„). For the two-component case 

(2.35) +J2 Ur,0,, (D^Fi) + J2 V.P., (D'^^^^) 

^ OUn+k ^ OVii+k 

Example 2.21. The KvM lattice H2.3|) has recursion operator 
TZ = u„(l + D)(u„D-D-i«„)(D-l)-i— I 

Un 

(2.36) = u„D"^ + (7i„ + w„+i)I + u„D + u„(u„+i - u„_i)(D - I)"^ — I. 

Un 

Example 2.22. The Toda lattice H2.4(l has recursion operator 

-Unl -D-l - I + (Vn-l - Vn){B - 1)-^^ I 



(2.37) TZ 



-w„I - v„D -M„+iI + t;„(u„ - ■u„+i)(D - I) ^— 1 



In Section we will give an algorithm for the computation of scalar recursion 
operators like H2.36|l . In Section [S] we cover the matrix case and show how H2.37|l is 
computed. The algorithms complement those for recursion operators of PDEs pre- 
sented in 'HGQQJ and elsewhere in these proceedings •BHSQ4! . We now introduce 
two powerful tools which will be used in the computation of densities and fluxes. 

3. The Discrete Variational Derivative (Euler Operator) 

Definition 3.1. A function i?„(u„) is a total difference if there exists another 
function J„(u„), such that En = AJn = (D — I) J„. 

Theorem 3.2. A necessary and sufficient condition for a function En with 
positive shifts up to level po, to be a total difference is that 

(3.1) C^^UEn) = 0, 
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where /^u°2 the discrete variational derivative (Euler operator) |ASYOd] defined 
by 

p. PO p. 

(3.2) = 7r-£ ^'^^ = -^(I + + + . . . + D-P°). 

k—u 

A proof is given in e.g. jHH03| . 

Remark 3.3. To verify that an expression E(un-q, ■ ■ ■ ,Un, - ■ ■ , Un+p) involving 
negative shifts is a total difference, one must first remove the negative shifts by 
replacing En by E„ — D^-Bn- Applied to En, H3.2|l terminates a,t po — p + q. 

We now introduce a tool to invert the total difference operator A = D — I. 

4. The Discrete Homotopy Operator 

Given is an expression En (free of negative shifts). We assume that one has 
verified that En € KerC^}, i.e. £ui(i?„) — 0. Consequently, En S ImA. So, E1J„ 
such that En — AJn- To compute J„ — A~^(i?„) one must invert the operator 
A = D — 1. Working with the formal inverse, 

(4.1) A"i =D"i+D-2 + D-3 + ..., 

is impractical, perhaps impossible. We therefore present the (discrete) homotopy 
operator which circumvents the above infinite series. In analogy to the continuous 
case |093j . we first introduce the discrete higher Euler operators. 



Definition 4.1. The discrete higher Euler operators are defined by 

d 



k—i 



The higher Euler operators all terminate at some maximal shifts pi. 
Example 4.2. For scalar component it„, the first higher Euler operators are: 

(4.3) 4°) = -^(I + + + + ■ ■ • + B-P«), 

OUn 

(4.4) 4) = -^(D-i+2D-2 + 3D-3 + 4D-4 + ...+p,D-fi), 

OUn 

(4.5) 42) = ^(D-2 + 3D-3 + 6D-4 + i0D-5 + --- + ip2(p2-l)D-''^). 

OUn 2 

Note that (|4.3|) coincides with (|3.2|) . Similar formulae hold for cl^} . 
Next, we introduce the homotopy operators. For notational simplicity we show 
the formulae for the two component case u„ = (u„,i;„). 



(4.6) / (ii,„(u„)[Au 

n\ J2,n 

(u„)[Au„])- 

Jo 



Definition 4.3. The total homotopy operator is defined as 

dX 
T' 

where the homotopy operators are given by 

Pi-i P2-1 
(4.7) ji,„(u„) = ^ (D - I)'(ii„4+i)), j2,n(u„) = E (D - I)^K4+'^)- 

i=0 i=0 
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Theorem 4.4. J„ = A^^iEn) can be computed as J„ — Ti.{En). 

A similar theorem and proof for continuous homotopy operators is given in 
|093| . By constructing a similar variational bicomplex the theorem still holds in 
the discrete case. See |HM02L IMH02| and elsewhere in these proceedings |M04| . 

Remark 4.5. Since the theory for the recursion operator can be formulated 
in the language of Lie algebra complexes, the results in |SW01a] are immediately 
applicable. This gives one explicit conditions from which the existence of infinite 
hierarchies of local symmetries and cosymmetries can be concluded. Also, the 
problems with the definition of a recursion operator as treated in ISWOlbj also 
play a role in the discrete case. 

Remark 4.6. Note that pi and p2 in H4.7() are the highest shifts for u„ and 
in En- Furthermore, jr,n(u„)[Au„] means that in jr,n(u„) one replaces u„ Au„, 
and, of course, u„+i Au„+i, u„+2 ^ Au„+2, etc. 

Remark 4.7. In practice, one does not compute the definite integral 1)4. 6|l since 
the evaluation at boundary A = may cause problems. Instead, one computes the 
indefinite integral and evaluates the result at A = 1. 



5. Algorithm to Compute Densities and Fluxes 

As an example, we compute the density-flux pair H2.21|l for (|2.4|l . Assuming 
that the weights (|2.8I) are computed and the rank of the density is selected (i? = 3 
here), the algorithm has three steps. 
Step 1: Construct the form of the density 

List all monomials in u„ and u„ of rank 3 or less: ^ = {un'^, u„^, u„w„, u„, f„}. 

Next, for each monomial in introduce the correct number of t-derivatives so 
that each term has weight 3. Thus, using 1)2. 4|l . 

d\ 3^ 3 

d 2\ r, • r> r, ^ , , . 

— (U„ j 2UnUn = 2u„t;„_i - 2UnVn, " ^ ^ """"" ~ "n+l^^m 

d'/ ^ d ^ . ^ d , 

Gather all terms in H = {uj^ ,UnVn-i,UnVn,Un-iVn-i,Un^iVn}. Identify members 
belonging to the same equivalence classes and replace them by their main represen- 
tatives. For example, UnVn-i = Un+i^n, so both are replaced by UnVn-i- So, Ti is 
replaced by Z = {u„^, u„w„_i, u„w„}, which has the building blocks of the density. 
Linear combination of the monomials in I with constant coefficients Ci gives 

(5.1) Pn ^ ClUn^ + C2UnVn-l + C^UnVn- 

Step 2: Determine the coefRcients 

Require that (|2.13|) holds. Compute Dtp„. Use (12. 4|) to remove m„ and w„ and their 
shifts. Thus, 

En = Dfp„ = (3ci - C2)w^W„_i -I- (C3 - icijU^^Vn + (c3 - C2)Vn-lVn 

(5.2) +C2Un-lUnVn-l + C2Wn_i - C3U„U„+lW„ - CaU^. 
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Compute En ~ Di^n. First, apply H3.2|l for component u„ to En ■ 

OUn 

= 2(3ci - C2)u„w„_i + 2(c3 - 3ci)m„w„ 

(5.3) +(C2 - C3)u„_ii;„_i + (C2 - C3)Un+lVn. 

Second, apply (|3.2|l for component Vn to _E„ : 

= (3ci - C2)u^j+i + (C3 - C2)w„+1 + (C2 - C3)u„M„+i 

(5.4) +2(c2 - C3)w„ + (c3 - 3ci)m^ + (03 - C2)w„_i. 
Both H5.3|l and H5.4|l must vanish identically. Solve the linear system 

(5.5) S = {3ci - C2 = 0, C3 - 3ci = 0, C2 - C3 = 0}. 

1, into H5.1|l 

1 



The solution is 3ci ~ C2 ~ c^. Substituting 01 — ^,02—0-3 



(5.6) 



-Un^ + Un{Vn-l + V^) ■ 



Step 3: Compute the flux 

In view of (|2.13|) . we will compute Jn = A~^(—En) with the homotopy operator 
introduced in Sectional Alternative methods are described in (GH98. .HH03) . 
Insert ci = ^,02 = C3 — 1 into (|5.2|l and compute 

(5.7) En = ^En = UnUn+lVn + - Un+lUn+2Vn+l - W^+i. 

We apply formulae (|4.7|l to —En- The pieces are listed in Tables ^ and El 



i 


c\::'\-En) 


{D-l)HunC\:+'\-En)) 



1 


Un-lVn-l+Un+lVn 


UnUn-lVn-l+UnUn+lVn 
Un+lUnVn-UnUn-lVn-1 



Table 1. Computation of ji,n(u„)(—£'„) 



i 


ci:,:'\~En) 


{B-iy{vnC[:+'\-En)) 








Table 2. Computation of j2,n(u„)(—£'„) 



Adding the terms in the right columns in Tableland Table [3 

(5.8) jl,n(u„)(-i?„) = 2w„U„+iU„, j2,«(u„)(-^„) = UnUn+lVn + 2vl_. 

Thus, the homotopy operator (|4.6|l gives 



Jn 



(5.9) 



(ji,n(u„)(--E„)[Au„] + j2,«(u„)(-£'„)[Au„]) ^ 



(3 A UnUn+lVn + 2Au„) d\ 
UnUn+lVn + W^. 
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After a backward shift, J„ = D ^ J„, we obtain the final result: 

(5.10) p„ = i + U„(u„_i + U„), Jn = Un-lUnVn-l + vl_j^. 

6. Algorithm to Compute Generalized Symmetries 

As an example, we compute the symmetry H2.30|l of rank (3,4) for (|2.4(l . The 
two steps of the algorithm ' GH99| are similar to those in Section 
Step 1: Construct the form of the symmetry 

Start by listing all monomials in u„ and v„ of ranks 3 and 4, or less: 

As in Step 1 in Section|3 for each monomial in K-i and /C27 introduce the necessary 
t-derivatives so that each term has rank 3 and 4, respectively. At the same time, 
use H2.4|l to remove all t— derivatives. Doing so, based on ICi we obtain 

(6.1) £1 = {ul,Un-lVn-l,UnVn-l,UnVn,Un+lV„}. 

Similarly, based on IC2, we get 

^2 = {Un,u'^_iVn-l,Un~lUnVn^i,u1vn_i,Vn-2Vn-l,v'^_j^,u1vn, 

(6.2) UnUn+lVn,ul^iVn,Vn-lVn,vl,V„Vn+l}. 

In contrast to the strategy for densities, we do not introduce the main represen- 
tatives, but linearly combine the monomials in Ci and C2 to get the form of the 
symmetry: 

Gi = Ciul^ + C2Un-lVn-l + CsUnVn-l + CiUnVn + C5Un+lVn, 
G2 = Ceu'^ + CrU^^iVn^i + C8Un~lUnVn-l + Cgul^Vn-l + CioVn-2Vn-l 
+ C11 vl_i + C12 ul^Vn + C13 UnUn+lVn + C14 U^n+l^n + Cl5 Vn-lVn 

(6.3) -|-Ci6W^ + Ci7W„D„+i, 

with constant coefficients c^. 

Step 2: Determine the coefficients 

To compute the coefficients Ci we require that l|2.22|l holds on solutions of 1)2. 1|1 . 
Compute DjCi and DtG2 and use l|2.4|l to remove Un,Vn, and their shifts. This 
gives the left hand sides of H2.22|l . 

Use (j^^ to get the right hand sides of (j??^ : 

i^i'(u„)[G] = D-IG2-IG2, 

(6.4) F2'(U„)[G] = VnlGi ~ VnDGi + {Un - Un+l)lG2. 

Substitute (|6.3() into (|6.4|l and equate the corresponding left and right hand sides. 
Since all monomials in Un,Vn and their shifts are independent, one obtains the 
linear system that determines the coefficients The solution is 

Cl = Cg = C7 = C8 Cg = Cio = Cii = C13 = C16 = 0, 

(6.5) C2 = C3 = -C4 = -C5 = C12 = -Ci4 = Ci5 = -C17. 

With the choice C17 = 1, the symmetry (I6.3|l finally becomes 

Gl = VniUn + Un+l) - Vn-l{Un-l - Un), 

(6.6) G2 = Vn{u1^^ - u'^ + Vn+i - Vn-l). 
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7. Algorithm to Compute Scalar Recursion Operators 

In this section we show how to compute the recursion operator H2.36|l of H2.4|l . 

Again, we will use the concept of rank invariance to construct a candidate 
recursion operator. The defining equation H2.32|l is then used to determine the 
coefficients. 

We observe that (|2.36|l in expanded form naturally splits into two pieces: 

(7.1) 7^ = 7^o + 7^l, 

where TZq contains only terms with shift operators D~^, I, and D, and TZi has terms 
involving (D — 

Step 1: Determine the rank of the recursion operator 

In view of H2.31|l and assuming that the symmetries are linked consecutively {s — 1), 

the recursion operator TZ has rank 

(7.2) R = rank 7^ = rank G^^^ - rank G^^^ =3-2 = 1, 

where we used ()2.27(l and (|2.28(l . If the assumption turns out to be correct, then 
the recursion operator has rank 1. If the assumption were falls because symmetries 
are not linked consecutively, then R must be adjusted and the subsequent steps 
must be repeated. See jBHS04] for examples of PDEs for which that happens. 
Step 2: Determine the form of the recursion operator 
We split this into two sub-steps. 

(i) Determine the form of TZq 

The candidate TZq is a sum of terms involving D~^,I, and D. The coefficients of 
these terms are linear combinations of Un-i,Um and Un+i, so that all terms have 
the correct rank. 

TZo = (C1U„-1 + C2U„ + C3U„+l) D"-^ + (C4U„_1 + C5U„ + C6U„+l) I 

(7.3) -|-(C7U„_1 + CsUn + C9U„+i) D, 

where the Ci are constant coefficients. 

A few remarks are in place. First, in TZq we moved the operators D^^,I, and 
D all the way to the right. Second, the maximum up-shift and down-shift operator 
that should be included can be determined by comparing two consecutive symme- 
tries. Indeed, if the maximum up-shift in the first symmetry is Un+p, and the max- 
imum up-shift in the next symmetry is u,i+p+ri then TZq must have D, D^, • ■ • , D*". 
The same line of reasoning determines the minimum down-shift operator to be in- 
cluded. In our example, there is no need to include terms in D~^,D^, etc. Third, 
the coefficients of the operators can be restricted to linear combinations of the terms 
appearing in F. Hence, no terms in Un±2, Un±3 and so on occur in 17.3|l . 

(ii) Determine the form of TZi 

As in the continuous case |HG99j , TZi is a linear combination (with constant coef- 
ficients Cjk) of sums of all suitable products of symmetries and covariants (Frechet 
derivatives of densities) sandwiching (D — I)~^, i.e. 

(7.4) 7^l=5^5^5,,G(^■)(D-I)-Vi'=''. 

J k 

For H2.3|l . G^^^ in H2.27|l and = ln(M„) in H2.15|l are the only suitable pair. 
Indeed, using H2.24|l we have pn^ = ln(u„)' = ^I, which has rank -1. Combined 
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with G*-^-* of rank 2, we have a term of rank 1. Other combinations of symmetries 
and covariants would exceed rank 1. Therefore, 

(7.5) Til = ciom„(m„+i - u„_i)(D - I)"^(— )I 
where cio a constant coefficient. Using (|7.1|l and renaming cio to cio, 

7^ = (ciU„_l + C2U„ + C3U„+l)D"^ + (C4U„_1 + C5U„ + C6U„+l)I 

(7.6) +(C7U„_1 + C^Un + C9U„+i)D + CioUn{Un+l " M„-l)(D - I)"^( — )I 

Un 

is the candidate recursion operator for (|2.3|l . 
Step 3: Determine the coefficients 

Starting from (|7.6|l . we use H2.34|l with F = u„(m„+i — u„_i) to compute 7?.'[F]. 
The partial results will not be shown due to length. 
Using (|2.23(l . we compute 

(7.7) F' = -U„D-1 + [Un+l - Un-l)l + u„D. 

Then we compute Tio F' and i^' o Ti. After substituting the pieces into (|2.82|l we 
simplify the resulting expression using rules such as 

(D-I)"^D = D(D-I)"^ =I+(D-I)^\ 

(7.8) (D-I)-iD-i = D-i(D -I)-i = -D-i + (D -I)-i. 
We further simplify by recursively using formulas like 

D[/(U„)(D-I)-V(W„)I = U{Un+l)V{Un)l+U{Un+l){D-l)-^V{Un% 

(7.9) {B~l)-^U{Un)V{Un)B - U{Un-l)V{Un-l)l+{B-l)-^U{Un-l)V{Un-l)l. 

Finally, we equate like terms to obtain a linear system for the q. Substituting the 
solution 

(7.10) Ci C3 = C4 = C7 = Cg = 0, C2 = C5 = C6 = C8 = CiQ = 1 

into (|721) we obtain the final result 

(7.11) TZ = U„D"^ + {Un + M„+l)I + U,iD + Un[Un+l - U„-l)(D - 1)^^—1. 

Un 

A straightforward computation confirms that TJ-G^^-* = G*-^^ with G^^^ in H2.27|l and 
G(2) in ^T^. 

8. Algorithm to Compute Matrix Recursion Operators 

We construct the recursion operator H2.37|l for (|2.4|l . Now all the terms in (|2.32|l 
are 2x2 matrix operators. 

Step 1: Determine the rank of the recursion operator 

The difference in rank of symmetries is again used to compute the rank of the 
elements of the recursion operator. Using l|2.8(l . 1)2.29(1 and ((2.30(1 . 

(8.1) rankG(i) = ( 3 ) ' I'^i^^G^^) = 

Assuming that TiG'^^^ = G^^), we use the formula 

(8.2) rank 7^y = rank Gf - rank cf^ , 
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to compute a rank matrix associated to the operator 

1 



(8.3) rank7^ = 



2 1 



Step 2: Determine the form of the recursion operator 

As in the scalar case, we build a candidate TZq : 

' (7eo)ii (7^o)l2 



(8.4) 7^o = 

with 



{no)2i (7^o) 



22 



{no)i2 = C3D-1 + C4I, 

{1^0)21 = (C5U^ + CeUnUn+l + CjU^n+l + Cs^n-l + CgWn) I 

+ (cioM^ + ciiu„u„+i + ci2U^+i + Cl3^;„_l + ciaVu) D, 

(7?-o)22 = (cisUn + Ci6W„+l) I. 

Analogous to the scalar case, the elements of matrix TZi are linear combinations 
with constant coefficients of all suitable products of symmetries and covariants 
sandwiching (D — I)~^. Hence, 

(8.5) ^^5,feG«(D-l)-i®p«', 

j k 

where (8) denotes the matrix outer product, defined as 



G« i ' ^y'-'' '^"''^"l g(^")(d-i)-vS' g(^")(d-i)-V^^^' 



n.2 



Only the pair {G'^^\ pl°^' ) can be used, otherwise the ranks in (|8.3|l would be ex- 
ceeded. Using H2.26|l and (|2.f8(l we compute 

(8.6) pif-{0 1^:1), 

Therefore, using (|8.5|l and renaming cio to cij, 



(8.7) 7^l = 



ci7(i;„-i -i;„)(D-I)-i-^l 



CuVn {Un - Un+l) (D - I) ^ ;^ 1 

Adding (|8.4|l and (|8.7|l we obtain 

■ 7^n 7^l2 

T^21 T^22 



with 



7^11 = (C1U„ + C2U„+l)l 

7^l2 = C3D-i+C4l + ci7(i^„-i -i^„)(D-I)-i— 1, 

7^21 = (C5U^ + C6W„M„+l + CyU^j+i + C^Vn-l + CgWn) I 

-|-(cioU^ + Cii-U„-U„+i + Ci2M^+i + Ci3t;„_i + CiiVn) D, 

7?-22 = (C15U„ + Ci6U,i+l)r+ ci7i;„(u„ - U„+i)(D - I)"^ — f. 
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Step 3: Determine the coefficients 

All the terms in H2.32|l need to be computed. The strategy is similar to the scalar 
case, yet the computations are much more cumbersome. Omitting the details, the 
result is: C2 = C5 = cg = cy = Cg = cio = Cn = C12 = C13 = C15 = 0, ci = C3 = 
C4 = cg = ci4 — C16 = — 1, and C17 — 1. Substitution the constants into H8.8|) gives 

(8 9) n=( """^ -D-i-I+K_i-i;„)(D-I)-i^I 

V -«rJ-w„D -w„+iI + i;„(m„ -u„+i)(D-I)"i 

It is straightforward to verify that TZG^'^^ = G^^) with G^^'> in itT^ and G'^) in 

9. More Examples 
Example 9.1. The modified Volterra lattice |ASYOOL IHHOS] , 

(9.1) iln = M^i(m,i+1 - Un-l), 

has two non-polynomial densities pn^ = ^^j- and pn"^ ~ ln('u„), and infinitely many 
polynomial densities. The first two symmetries, 

(9.2) G« = ul{ur,+i-u„^i), 

(9.3) G^^) = ulul^^{Un+Un+2) -ul_-^ul{Un-2+Un), 

are linked by the recursion operator 

(9.4) TZ = w^D-i + 2ii„u„+il + ulY) + 2ul{u„+i ~ w„_i)(D - 1)-^— I. 

Un 

Example 9.2. The AL lattice (|2.2|) has infinitely many densities GH98 and 
symmetries ^H99 . The recursion operator is of the form H8.8|l with 

Un = F„D-i -u„A-i«„+iI-u„_iF„A-i^I, 

TZi2 = -u„u,i-il - u„A~^u„_iI - u„_iP„A~^^I, 

7^21 = w„i'„+iI + w„A~^w„+iI + w„+iP„A"^^I, 

(9.5) 7?.22 = (m„w„+i + u„_iu„)I + P„D + u„A"^u„_iI + Vn+iP^A"^^!, 

where P„ = 1 + m„u„ and A = D — I. 

This recursion operator has an inverse, which is quite exceptional. The inves- 
tigation of the recursion operator structure of the AL lattice is work in progress. 
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